Mean-field phase diagram for Bose-Hubbard Hamiltonians with random hopping 
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The zero-temperature phase diagram for ultracold Bosons in a random ID potential is obtained 
through a site- decoupling mean-field scheme performed over a Bose-Hubbard (BH) Hamiltonian 
whose hopping term is considered as a random variable. As for the model with random on-site 
potential, the presence of disorder leads to the appearance of a Bose-glass phase. The different 
phases -i.e. Mott insulator, superfiuid, Bose-glass- are characterized in terms of condensate fraction 
and superfiuid fraction. Furthermore, the boundary of the Mott lobes are related to an off-diagonal 
Anderson model featuring the same disorder distribution as the original BH Hamiltonian. 

PACS numbers: 03.75.Lm, 05.30.Jp, 64.60.Cn 



I. INTRODUCTION 

In the framework of ultracold atom physics, the exper- 
imental tunability of control parameters pertaining each 
model Hamiltonian has provided a powerful tool to in- 
vestigate situations of fundamental physical interest [l| . 

One of the most intriguing features about ultracold 
atoms is the possibility to engineer a defect-free periodic 
potential, as opposed, for instance, to the typical frame- 
work of solid-state physics. However, on the one hand the 
interplay between disorder and interactions in Bosonic 
systems has attracted much theoretical attention since 
the seminal work by Fisher [1| , and on the other several 
techniques such as laser speckle field [3j], the superpo- 
sition of different optical lattices with incommensurate 
lattice constants have proven the experimental 

relevance of disordered systems of ultracold atoms. 

In the present paper we will deal with the effect of dis- 
order on the zero-temperature phase diagram of bosonic 
atoms loaded onto a ID lattice whose properties can be 
described in terms of Bose-Hubbard Hamiltonian. In par- 
ticular, we will focus on the case where the disorder af- 
fects the hopping term 



M 



^ ^ J 711 .in' h-C, 



(1) 



where a^^ {dm) represents the creation (destruction) 
operator on site m. The Hamiltonian parameters U, 
Jm.m' represent the two-body interaction and the (ran- 
dom) hopping amplitude between neighboring sites re- 
spectively. 



In the recent past, many authors have approached the 
analysis of the disordered BH model with various tech- 
niques such as field-theoretic approaches [1, 0, H, [sl, 
decou pling (or Gunzwiller) mean-field approximations 
^ [lOl. [Til [l3 . [Tsj . quantum Monte-Carlo simulations 
IllQ, among many others, e.g. [H, [13, [H, [iSj, l20| . 

Following (2]| , we employ a site-decoupling mean-field 
approximation (SDMFA), which allows to capture all of 
the essential features of the phase diagram of the model 
([1]). The phases of the zero-temperature phase diagram 
are determined through the calculation of two differ- 
ent observables: the condensate fraction, defined as the 
largest eigenvalue of the one-body density matrix and the 
superfiuid fraction, defined as the system response to the 
coupling to an external field [1, H^l ■ 

At zero temperature, as already pointed out in for 
the on-site disordered BH model, we expect the presence 
of three phases: the Mott-Insulating (MI) phase, where 
both condensate fraction and superfiuid fraction are zero, 
the superfiuid phase, where both superfiuid phase and 
condensate fraction are different from zero, and, finally, 
the Bose-Glass phase, which has zero superfiuid fraction 
and finite condensate fraction, which represents a typical 
feature of disordered hopping systems. 

In Section |TT1 we introduce the site decoupling mean- 
field approximation for the case given by Hamiltonian ([T]) 
and we depict the zero-temperature phase diagram, dis- 
cussing similarities and differences between the random- 
hopping and the random on-site potential case. 
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In Section IIIIl we discuss the stability of the Mott 
phase through the stability analysis of the recurrence 
map induced by the SDMFA. This analysis will be per- 
formed comparing the results obtained by numerical 
exact diagonalization and analytical results based on 
random-matrix theory. 
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II. SITE-DECOUPLING MEAN-FIELD SCHEME 

The SDMFA was introduced in Ref. [l0[ , and relies on 
the approximation 

alnam' = al^'ttm + flmam' " Ctmam' (2) 

where the am 's are mean- field variables to be determined 
self-consistently. The above posit turns the BH Hamil- 
tonian ^ into a mean-field Hamiltonian that is the sum 
of on-site contributions. 



n = 5^7^, 



(3) 



(4) 



where rim = oL'^m the usual bosonic number operator 
and the disorder related to the hopping term has been 
embedded into the adjacency matrix Am,m' ■ For nearest- 
neighbor hopping on a ID system the adjacency matrix 
can be written as 



',— l^m'.m—l ~1~ 6 ^ Sm^m' ,m+l (5) 



where (j) takes into account a possible coupling to an ex- 
ternal field (Peierls factors), Sm G ^ accounts for a pos- 
sible inhomogeneity of the hopping amplitude. Here we 
consider a ID lattice comprising TV sites with periodic 
boundary conditions, i.e. 



sn+1 — Si, So 



sn- 



(6) 



The ground state of Hamiltonian ^ is clearly a prod- 
uct of on-site states, 



I*) 



(7) 



and the requirement that the relevant energy is mini- 
mized results in the self-consistency constraint [l^ 



rn 771' *~*-m i 



Ctm = (V'mlamlV'm) (8) 



The decoupling mean-field approach has proved to pro- 
vide satisfactory qualitative phase diagrams for homo- 
geneous lattices [13, [H, [l^] and superlattices 

Mm. 



A. Numerical Simulation 

In the present section we report the results of the appli- 
cation of SDMFA to the zero-temperature phase diagram 
calculation. 
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FIG. 1: Values of taken into account for the determi- 
nation of the phase diagrams. Upper plot has weaker dis- 
order amplitude (Amax = 0.1) compared to the lower one 
(A^ax = 0.2). 



For our calculations we have considered a lattice 
formed by 100 sites, with values of the adjacency ma- 
trix given by 

+ (1 + '^m+l )Sm' ,m+l (9) 

where Am is an uncorrelated random variable uniformly 
distributed between Amax and —Amax, with the con- 
straint Aniax < 1 in order to preserve Jm,m' > 0. In 
Fig. Ill Al we have represented two explicit realizations 
of disorder taken into account for subsequent mean-field 
calculations of phase diagrams. 

The value of Amax has been kept small enough to en- 
sure reasonable self-averaging for the system under in- 
vestigation. Higher disorder amplitudes would require 
larger chains in order to obtain results independent of 
the specific disorder realization. 

The different phases have been identified through the 
evaluation of two observables, namely the condensate 
fraction, defined as the largest eigenvalue of the one-body 
density matrix 



fc — ^^-^{Pm-m') 1 Pm.m' — (^rr 

and the superfluid fraction defined as 

E{^) - EiO) 



fs 



lim 

0^0 



J{N)4>'^ 



^)/N, (10) 



(11) 



where E((j) is the ground-state energy corresponding to 
the presence of a Peierls phase (j) in the hopping term 

The MI phase is characterized hy fc = fs — 0, the 
superfluid phase (SF) by ^ and fs ^ 0, while the 
phase where /c ^ and fs — is recognized as the Bose- 
glass (BG) phase 2]. 

In absence of disorder, the variation of the control pa- 
rameters - chemical potential and hopping amplitude- 
always leads to a transition from a phase with both fs 
and fc equal to zero to a phase with fs and fc different 
from zero, excluding then the presence of a Bose-glass 
phase. 

On the other hand, if on-site disorder is present, the 
MI phase is (possibly) separated from the SF phase by 
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a BG phase. Likewise, when disorder affects the hop- 
ping term alone a BG crops up, as it is shown in Fig. 
Ill Al However the distribution of the BG phase in the 
parameter space is quahtatively different. For example, 
with on-site disorder MI and SF phase are separated by 
a BG phase as J goes to zero, while with disorder on the 
ho ppi ng term the BG phase tends to disappear for small 
J [2l|. This region of the phase diagram seems to be 
a good starting point for future investigations by means 
of a cluster MF approach [2^, [13] , because it would re- 
veal possible MI phases that are not detectable through 
the single-site mean-filed technique implemented in the 
present paper. 

In Fig. Ill Al we have represented the first lobe of the 
zero temperature phase diagram as obtained by SDMFA 
for Amax = 0.1 and A^ax = 0.2. 



III. STABILITY OF THE MOTT PHASE 

In the present Section we introduce a procedure to de- 
termine the border of the MI phase which is based on 
the stability of the self-consistency map induced by the 
mean-field procedure. 

As a first consideration, it is possible to state that the 
condition 7/c = for every k corresponds to the gapped 
insulating phase of the mean-field Hamiltonian[31 Indeed, 
in this case the local ground-states ([7]) are eigenvectors 
of the local number operator a^a 



n ^ — , 



(12) 



where [x] denotes the smallest integer larger than x, 
and hence the mean-field ground-state (O is a pure Fock 
state. As for the ordered Bose- Hubbard Hamiltonian, the 
relevant on-site energy is 



U 



n{n — 1) + /zn 



(13) 



Hence {i}}\a\il)) = at every site, and the sclf-consistcncy 
constraint ^ is satisfied. In other terms 7fe = is a fixed 
point of the map defined by Eq. ([8]). The gapped insu- 
lating phase is stable as long as this fixed point is stable. 
This is true as long as the magnitude of the maximal 
eigenvalue of the matrix A appearing in the linearized 
version of Eq.® is smaller than 1, see 114)) . In order 
to determine A, we assume |7fc| <C 1 and consider the 
(mean-field) kinetic term in Hamiltonian ^ as pertur- 
bative. If first order perturbation theory is performed 
one gets 



where 



F{x) 



(W -x){x- \x\) 



(14) 



(15) 
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FIG. 2: Plot of the zero-temperature phase diagram for 
Amax = 0.1 (upper panel) and Amax = 0.2 (lower panel). In 
black we have plotted the border between the MI and the 
BG phase, while the colored background corresponds to the 
superfluid fraction (grey=0). It is possible to notice how the 
BG phase covers an increasing area as Amax is increased. 



Hence the linearized version of the self-consistency map 
Eq. ([5]) can be written as 



J 



(16) 



where the matrix A is proportional to the adjacency ma- 
trix A: 



A = F U. 



(17) 
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Recalling the criteria for the stability of linear maps, the 
fixed point (a,„} — (equivalent to 7^ = 0) is stable 
whenever 



J 1 



(18) 



where Amax the eigenvalue of A with the largest magni- 
tude. 

The problem of determining the maximal eigenvalue 
of the matrix A can be reduced to the calculation of the 
maximal eigenvalue of the ( tridiagonal ) adjacency ma- 
trix A. Note that A is basically the one-particle Hamil- 
tonian for a non-interacting off-diagonal Anderson model 
(random hopping model). If Amax is the maximal eigen- 
value of the adjacency matrix, Eq. (fT5)) can be recast 
as 



J _ 1 

Tf-Jx 



(19) 



We have dealt with the calculation of the eigenvalues 
of A both numerically and analytically. The analytical 
approach consists in the determination of the spectral 
density of the matrix, given the probability distribution 
of the elements of the matrix following the approach out- 
lined by Dyson for and harmonic-oscillator chain fsij ]. 
while the former simply consists in the direct numerical 
evaluation of the matrix eigenvalues. 



A. Numerical analysis 

In Fig. [3] we report the maximum eigenvalue Amax 
for the random adjacency matrix A as a function of the 
strength of disorder for both a single realization and a 
disorder average over 1000 samples. The two panels re- 
fer to different lattice sizes. In every case, it is evident 
that the maximal eigenvalue Amax grows with increasing 
disorder strength, Amax • 

The stability condition given by Eq. (fT51) . along with 
the above considerations about the Amax-dependence 
from Amax are in agreement with the considerations of 
Section lll Al . as far as the MI phase is concerned. In par- 
ticular Amax can be thought of as a shrinking factor for 
the MI lobe when compared to a homogeneous situation. 



B. Analytical solution 

In this section we would like to describe an analytical 
method to obtain the largest eigenvalue of a (possibly in- 
finite) random adjacency matrix whose entries are given 
by Eq. ([9]). For some particular disorder distributions 
it is possible to carry through the analytical calculation 
and, as a consequence, solve Eq. (I19p without finite-size 
effects. 
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FIG. 3: Amax-dependence of the largest eigenvalue of the ad- 
jacency matrix A, for 200 sites (upper panel) and 100 sites 
(lower panel). By comparison between the single realization 
plots (full lines) and the averaged ones (dashed lines), it is 
possible to see how, increasing the lattice size, the effect of 
the specific randomness realization is decreased 
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FIG. 4: Boundary of the MI region for different values of 
A 

max . 



It is worth mentioning that, in principle, the informa- 
tion gained through this approach is richer. In fact for 
some specific realizations of disorder it is possible to ob- 
tain the full integrated density of states M{z) and not 
simply the largest eigenvalue of the adjacency matrix. 

The solution of the problem can be obtained following 
the approach proposed by Dyson for the solution of a 
linear chain of harmonic oscillators. 

Here we will sketch Dyson's approach, outlining the 
connection with our problem, providing the portion of 
M{z) needed to obtain the largest eigenvalue in view of 

In^j^. the problem of a linear harmonic chain with 
springs of random elastic is reformulated as a tridiagonal 
matrix diagonalization problem. The matrix A has the 
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following form 



A 



J+ij 



-A,7.j+i = i^j 



1/2 



(20) 



which is related to the matrix A defined in Eq. ([5|) by a 
unitary transformation U{0) 



with 



and 



Uj,Ki9)=6j,KeM^dJ) 



A = u{e)Au{e) 



(21) 



(22) 



1 /2 

and setting sj — Xj in Eq. ([5]). The unitarity of U{9) 
allows to state that the eigenvalues of A are equal to the 
eigenvalues of A, hence the procedure followed in [sij can 
be directly mapped onto our problem. 

The core of this approach resides in the definition of 
the characteristic function of the chain 



r2(x) = lim 



N 



(23) 



where N is the size of the matrix, and ujj the desired 
eigenvalues of the matrix under consideration. 

The density of states D{z) and the integrated density 
of states 



M{z) = / D{z')dz' 



can be obtained from the characteristic function through 
the following relation 



D{z) = ^-Im 

having defined 



limO'(-- 

e^O Z 



(24) 



n'{x) = 



dx 



The determination of ^l{x) is obtained through a power 
series expansion 



oo 

n(x)= lim — V(-l)"-irr(A2"). 



(25) 



n=l 
I 2n\ 



The determination of Tr(A ") leads to the following re- 
lation 



1 ^ 



(26) 



a=l 



having defined ^(a) through the continued fraction 

e(a) =xAa/(l + xAa+i/(l + a;A,+2/(... . (27) 



If, as it may be safely assumed in our case, the vari- 
ous values of are independent random variables with 
probability distribution G(A), the variable 



xXa 



1 + S,ia + 1) 



will have probability distribution F(^) satisfying the fol- 
lowing integral equation 

m = r ^(^')^ + ^'/^)] (28) 

Jo ^ 



With the normalization condition 



we obtain 



/■oo 

/ = 1 

^0 

Jo 



If we assume a Poissonian form for G(A) 



G„(A) = 



D! 



A"-icxp(-nA), 



(29) 



(30) 



(31) 



Eq. has an analytical solution. Hence it is possible 
to obtain the integrated density of states in closed form 
in terms of integral functions. 

The integrated density of states for A can be simply 
obtained back by posing 



M'^iz) = M(z2) 



(32) 



since in [31[ M{z) is defined as the proportion of eigen- 
values for which Uj < z, while M^{z) as the proportion 
of eigenvalues with \ujj \ < z. As a simple check, we pro- 
vide here the expression for M^{z), corresponding to the 
case without disorder, i.e. n = oo in Eq. pip . 



M^{z) = iarccos [l - l/2z2] z<2 



1 



z> 2. 



(33) 



which, as expected, coincides with the well known result 
for a ID homogeneous system. 

On the other hand, if we are interested in the deter- 
mination of the maximal eigenvalue of A in presence of 
(weak) disorder whose distribution can be related to that 
expressed by Eq. ((3T|) . we can consider a large-n expan- 
sion of M^{z) for z > 2. The expression of in this 
case is given by 



M^{z) ~ 1 exp [—a — 2n (sinha — a)] 



(34) 



with a — arccosh (1/22:^ — l). Fig. [5] shows the behav- 
ior of the integrated density of states in the vicinity of 
the band edge at Amax, for three different values of the 
disorder strength. In Ref. [2l[ the behavior of the cor- 
responding density of states is related to the presence 
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FIG. 5: Comparison between three analytical solutions for 
M{Z) = 1, it is possible to see how, for decreasing disor- 
der (n increasing) the solution approaches the solution of the 
problem without randomness. In the x-axis the variable z 
has been normalized with respect to the maximum eigenvalue 

^max 
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FIG. 6: Plot of the "shrinking factor" for the MI lobe bound- 
ary as a function of the disorder intensity (1/n — > represents 
the case without disorder). 
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of a BG phase outside the MI region. In that case a 
possible direct transition from the MI to the SF phase 
is possible for small disorders and specific values of the 
chemical potential, and signaled by a singularity at the 
band edge in the density of states, similar to the Van 
Hove singularity characterizing the homogeneous case, 
Eq. (|33|) . Conversely, in the present case, where den- 
sity of states depends on the chemical potential through 
an overall multiplicative factor, an infinitesimal disorder 
is sufficient to smear the discontinuity in the density of 
states. Hence an intermediate BG phase is expected for 
every value of the chemical potential, which is in agree- 
ment with the previously noted shape of the BG phase 
in Figs. EAlindlll 

IV. CONCLUSIONS 



In this paper we have considered the effect of a ran- 
dom hopping term on the zero-temperature phase dia- 
gram of the Bose-Hubbard Hamiltonian. Analogously to 
what happens when a random on-site term is considered, 
we have observed the emergence of a Bose-glass phase. 
The analysis has been performed within a mean-field ap- 
proach both numerically and analytically. The bound- 
aries of the Mott lobes and the presence of a surrounding 
BG phase has been related to the spectral feature of an 
off-diagonal Anderson model. For the future, we plan to 
extend our research towards the finite-temperature case 
and to higher dimensions. 
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